#import the required libraries

!pip install srtm.py

!pip install geopandas;

!pip install pyshp;

!pip install rasterio;

!pip install keplergl

!pip install movingpandas

!pip install basemap

from descartes import PolygonPatch

import shapefile as shp

import matplotlib.pyplot as plt

import geopandas as gpd

import fiona

import numpy as np

import rasterio as rio

from rasterio.plot import show

import pandas as pd

import fiona.transform

import rasterio.sample

from scipy.spatial.distance import directed_hausdorff

from sklearn.cluster import DBSCAN

import movingpandas as mpd

from datetime import datetime, timedelta

import plotly.express as px

from mpl_toolkits.basemap import Basemap

from mpl_toolkits.mplot3d import Axes3D

import pylab

import math

import logging as mod_logging

import srtm

import srtm as mod_srtm

import urllib

import zipfile

import os

import scipy.io

import seaborn as sns

import matplotlib as mpl

import branca.colormap as cm

import folium

from scipy.stats import skewnorm

from folium import plugins

import geemap

%matplotlib inline

from pyvista import set_plot_theme

# Importing auxiliary libraries

import pandas as pd

import matplotlib.pyplot as plt

from IPython.display import Image

set_plot_theme('document')

mpl.style.use('default')

%matplotlib inline

import warnings

warnings.simplefilter("ignore")

# rayshader

library(rayshader)

library(tiff)

library(raster)

library(geoviz)

​

# chicago

options(warn=-1)

library(ggplot2)

library(tidyverse)

library(rgeos)

library(RColorBrewer)

library(maptools)

library(plotly)

library(scatterplot3d)

library(tidyverse)

library(RNHANES)

library(latticeExtra)

​

# spatial statistics

options(warn=-1)

library(spatstat)

library(raster)

%env HV_DOC_HTML=true

%load_ext rpy2.ipython

# reading the shapefile

fname = "//content//drive//MyDrive//GSP//academies-20160209-shp//academies-20160209.shp"

schools = gpd.read_file(fname)

schools.head()

fig, ax = plt.subplots(1, 1, figsize=(10,10))

schools.plot(cmap="OrRd",ax=ax)

plt.xlim(left=-5, right = 10)

plt.ylim(top=53, bottom = 40)

plt.title('French school districts')

plt.axis("off");

# importing data from geopandas datasets representing the earth

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))

world.crs # original projection

world.plot(figsize=(15, 10))

plt.title('World map from WGS84 (Latitude/Longitude)',fontdict={'fontsize':20})

plt.axis("off");

world_proj2 = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]

world_proj2 = world_proj2.to_crs("EPSG:3395") # world.to_crs(epsg=3395) would also work

mercator = world_proj2.plot(figsize=(15, 10));

mercator.set_title("Mercator Projection of The World Map",fontsize=20)

plt.axis("off");

#another testing of crs

#  check thos link for different crs: https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf

world_test = world.set_crs(32733, allow_override=True)

​

world_test = world_test.plot(figsize=(15, 10));

world_test.set_title("UTM (South) Projection of The World Map",fontsize=20)

plt.axis("off");

fpath = '//content//drive//MyDrive//GSP//srtm_germany_dsm//srtm_germany_dsm.tif'

​

def rasterio_open(f):

    return rio.open(f)

​

src_image = rasterio_open(fpath)

fig, ax = plt.subplots(1, figsize=(8, 8))

show(src_image, ax=ax,cmap='viridis')

plt.title('Germany raster map')

plt.axis("off")

plt.show()

file = "//content//drive//MyDrive//GSP//HYP_LR_SR_W//HYP_LR_SR_W.tif"

src_image = rasterio_open(file)

fig, ax = plt.subplots(1, figsize=(12, 12))

show(src_image, ax=ax,cmap='viridis')

plt.title('World raster map')

plt.axis("off")

plt.show()

dem = rio.open(file)

dem_array = dem.read(1).astype('float64')

​

dem_array = dem.read(1).astype('float64')

fig, ax = plt.subplots(1, figsize=(12, 12))

show(dem_array, cmap='Greys_r', ax=ax)

plt.title('World raster map')

plt.axis("off")

plt.show()

elevation_data = srtm.get_data(local_cache_dir="//content//drive//MyDrive//GSP//srtm")

print('CGN Airport elevation (meters):', elevation_data.get_elevation(50.8682, 7.1377))

image = elevation_data.get_image((500, 500), (45, 46), (13, 14), 300)

mod_logging.basicConfig(level=mod_logging.DEBUG,

                        format='%(asctime)s %(name)-12s %(levelname)-8s %(message)s')

​

geo_elevation_data = srtm.get_data()

​

image = geo_elevation_data.get_image((400, 400), (45, 46), (13, 14), 300)

image.show()

print('                Digital Elevation Model of Istra and Trieste')

geo_elevation_data.get_image((600, 600), (45, 46), (13, 14), 300)

# load file

str_name<-'H:/CY Tech Year 4/Semester 1/Geospatial/Datasets/finalmap.tif'

imported_raster=raster(str_name)

elmat = raster_to_matrix(imported_raster)

# create plot

elmat %>%

  sphere_shade(texture = "desert") %>% #add shade

  add_water(detect_water(elmat), color = "desert") %>% # add water

  add_shadow(ray_shade(elmat, zscale = 3), 0.5) %>% #add shadow

  add_shadow(ambient_shade(elmat), 0) %>% #add shadow

  plot_3d(elmat, zscale = 10, fov = 0, theta = 135, zoom = 1, phi = 45, windowsize = c(1000, 800))

Sys.sleep(0.2)

render_snapshot()

elmat %>%

  sphere_shade(texture = "desert") %>%

  add_water(detect_water(elmat), color = "desert") %>%

  add_shadow(ray_shade(elmat, zscale = 3), 0.5) %>%

  add_shadow(ambient_shade(elmat), 0) %>%

  plot_3d(elmat,

            solid = T,

            water = T,

            waterdepth = 0,

            wateralpha = 0.5,

            watercolor = "lightblue",

            waterlinecolor = "white",

            waterlinealpha = 0.5,

            fov=0,theta=135,zoom=0.75,phi=45, windowsize = c(1000,800))

Sys.sleep(0.2)

render_snapshot()

fname = './../Datasets/academies-20160209-shp'

french_schools = gpd.read_file(fname)

french_schools.head()

rng = np.random.default_rng()

french_schools["perc. exam success"] = rng.integers(low = 0, high = 100, size = len(french_schools))

x_map=french_schools.centroid.x.mean()

y_map=french_schools.centroid.y.mean()

print(x_map,y_map)

#checking color scheme 

colormap = cm.linear.YlGnBu_09.to_step(data=french_schools['perc. exam success'], method='quant', quantiles=[0,0.1,0.75,0.9,0.98,1])

colormap

map_ = folium.Map(location=[47, 3], zoom_start=5,tiles=None)

folium.TileLayer('OpenStreetMap',name="Light Map",control=False).add_to(map_)

colormap.caption = "% of exam success"

style_function = lambda x: {"weight":0.5, 

                            'color':'black',

                            'fillColor':colormap(x['properties']['perc. exam success']), 

                            'fillOpacity':0.75}

highlight_function = lambda x: {'fillColor': '#000000', 

                                'color':'#000000', 

                                'fillOpacity': 0.50, 

                                'weight': 0.1}

FRENCH_SCHOOLS=folium.features.GeoJson(

        french_schools,

        style_function=style_function,

        control=False,

        highlight_function=highlight_function,

        tooltip=folium.features.GeoJsonTooltip(fields=['name','perc. exam success'],

            aliases=['Name','% of exam success'],

            style=("background-color: white; color: #333333; font-family: arial; font-size: 12px; padding: 10px;"),

            sticky=True,zoom_start=0.1

        )

    )

colormap.add_to(map_)

map_.add_child(FRENCH_SCHOOLS)

map_

phd = pd.read_csv("./../Datasets/PhD.v3.csv")

phd.head()

phd_count = phd.groupby(["etablissement_rec"])["etablissement_rec"].count().reset_index(name = "count")

phd_count.head()

opendata = pd.read_csv("./../Datasets/opendata.csv",sep=";")[['uai - identifiant','Libellé','Géolocalisation']]

opendata.head()

opendata_phd_count = pd.merge(left=opendata, right=phd_count, left_on='Libellé', right_on='etablissement_rec')

opendata_phd_count.head()

CYlocation = list(map(float,str(opendata_phd_count[opendata_phd_count.Libellé=='CY Cergy Paris Université']['Géolocalisation'].values).split()[0][2:-2].split(',')))

#adding a marker on CY Cergy Paris Université location

folium.Marker(CYlocation, popup='CY Tech').add_to(map_)

#adding bubbles to the map based on the count of phds per institution

for i in range(0,len(opendata_phd_count)):

   coordinates = opendata_phd_count["Géolocalisation"][i]

   folium.Circle(

      location=coordinates.split(','),

      popup=opendata_phd_count.iloc[i]['Libellé'],

      radius=float(opendata_phd_count.iloc[i]['count']),

      color='blue',

      fill=True,

      fill_color='blue'

   ).add_to(map_)

map_

def kMedoids(D, k, tmax=100):

    # determine dimensions of distance matrix D

    m, n = D.shape

​

    np.fill_diagonal(D, math.inf)

​

    if k > n:

        raise Exception('too many medoids')

    # randomly initialize an array of k medoid indices

    M = np.arange(n)

    np.random.shuffle(M)

    M = np.sort(M[:k])

​

    # create a copy of the array of medoid indices

    Mnew = np.copy(M)

​

    # initialize a dictionary to represent clusters

    C = {}

    for t in range(tmax):

        # determine clusters, i. e. arrays of data indices

        J = np.argmin(D[:,M], axis=1)

​

        for kappa in range(k):

            C[kappa] = np.where(J==kappa)[0]

        # update cluster medoids

        for kappa in range(k):

            J = np.mean(D[np.ix_(C[kappa],C[kappa])],axis=1)

            j = np.argmin(J)

            Mnew[kappa] = C[kappa][j]

        np.sort(Mnew)

        # check for convergence

        if np.array_equal(M, Mnew):

            break

        M = np.copy(Mnew)

    else:

        # final update of cluster memberships

        J = np.argmin(D[:,M], axis=1)

        for kappa in range(k):

            C[kappa] = np.where(J==kappa)[0]

​

    np.fill_diagonal(D, 0)

​

    # return results

    return M, C

# from kmedoid import kMedoids # kMedoids code is adapted from https://github.com/letiantian/kmedoids

​

# Some visualization stuff, not so important

sns.set()

plt.rcParams['figure.figsize'] = (12, 12)

# Utility Functions

​

color_lst = plt.rcParams['axes.prop_cycle'].by_key()['color']

color_lst.extend(['firebrick', 'olive', 'indigo', 'khaki', 'teal', 'saddlebrown', 

                 'skyblue', 'coral', 'darkorange', 'lime', 'darkorchid', 'dimgray'])

​

def plot_cluster(traj_lst, cluster_lst):

    '''

    Plots given trajectories with a color that is specific for every trajectory's own cluster index.

    Outlier trajectories which are specified with -1 in `cluster_lst` are plotted dashed with black color

    '''

    cluster_count = np.max(cluster_lst) + 1

    for traj, cluster in zip(traj_lst, cluster_lst):

        if cluster == -1:

            # Means it it a noisy trajectory, paint it black

            plt.plot(traj[:, 0], traj[:, 1], c='k', linestyle='dashed')

        else:

            plt.plot(traj[:, 0], traj[:, 1], c=color_lst[cluster % len(color_lst)])

    plt.show()

dataset_link = '//content//drive//MyDrive//GSP//CVRR_dataset_trajectory_clustering//cross.mat'

data_folder = 'data'

filename = '%s/cross.mat' % data_folder

​

# Import dataset

traj_data = scipy.io.loadmat('//content//drive//MyDrive//GSP//CVRR_dataset_trajectory_clustering//cross.mat')['tracks']

​

traj_lst = []

for data_instance in traj_data:

    traj_lst.append(np.vstack(data_instance[0]).T)

# Plotting

​

for traj in traj_lst:

    plt.plot(traj[:, 0], traj[:, 1])

# 2 - Trajectory segmentation

​

degree_threshold = 5

​

for traj_index, traj in enumerate(traj_lst):

    hold_index_lst = []

    previous_azimuth= 1000

    for point_index, point in enumerate(traj[:-1]):

        next_point = traj[point_index + 1]

        diff_vector = next_point - point

        azimuth = (math.degrees(math.atan2(*diff_vector)) + 360) % 360

        if abs(azimuth - previous_azimuth) > degree_threshold:

            hold_index_lst.append(point_index)

            previous_azimuth = azimuth

    hold_index_lst.append(traj.shape[0] - 1) # Last point of trajectory is always added

    traj_lst[traj_index] = traj[hold_index_lst, :]

# 3 - Distance matrix

​

def hausdorff( u, v):

    d = max(directed_hausdorff(u, v)[0], directed_hausdorff(v, u)[0])

    return d

​

traj_count = len(traj_lst)

D = np.zeros((traj_count, traj_count))

​

# This may take a while

for i in range(traj_count):

    for j in range(i + 1, traj_count):

        distance = hausdorff(traj_lst[i], traj_lst[j])

        D[i, j] = distance

        D[j, i] = distance

​

k = 10 # The number of clusters

medoid_center_lst, cluster2index_lst = kMedoids(D, k)

​

cluster_lst = np.empty((traj_count,), dtype=int)

​

for cluster in cluster2index_lst:

    cluster_lst[cluster2index_lst[cluster]] = cluster

plot_cluster(traj_lst, cluster_lst)

mdl = DBSCAN(eps=400, min_samples=10)

cluster_lst = mdl.fit_predict(D)

​

plot_cluster(traj_lst, cluster_lst)

# import data

df_spacetime = pd.read_csv('//content//drive//MyDrive//GSP//MPIAB White Stork South African Population Argos.csv')

df_spacetime.head(2)

#initialise the figure object

fig = plt.figure(figsize=(14,10))

​

#initialise the Map using Basemap. You can manually set the area of the world you want to 

#be displayed. You can also get various variations of the maps.

​

m = Basemap(projection='ortho', lat_0=10, lon_0=13, resolution='i')

​

m.fillcontinents(color='#191919',lake_color='#000000') # dark grey land, black lakes

m.drawmapboundary(fill_color='#000000')                # black background

m.drawcountries(linewidth=0.1, color="w")              # draw countries

​

#extract the coordinates into the variables x and y

x, y = m(df_spacetime['location-long'].values, df_spacetime['location-lat'].values)

​

​

#plot the positions of the birds as they migrate using the scatter plot.

#Refer the documentation of Basemap for a better understanding of the various parameters

m.scatter(x, y, c="#1292db", lw=0, alpha=1, zorder=5, s = 1)

​

#title for the image

plt.title("Migration of Lesser Black-backed Gulls birds",fontsize=20)

​

#plot the image

plt.show()

set(df_spacetime['tag-local-identifier'])

def plot_3d(df):

    #initialise a figure object

    fig = plt.figure(figsize=(14,10))

    #add a 3D object

    ax = fig.add_subplot(111, projection='3d')

​

    groups = df.groupby('tag-local-identifier')

​

    # Plot

    # fig, ax = plt.subplots()

    # ax.margins(0.05) # Optional, just adds 5% padding to the autoscaling

    # for name, group in groups:

    #     ax.scatter(df['location-lat'], df['location-long'],df['event-id'], label=name, lw=0, alpha=1, zorder=5,c=df['tag-local-identifier'])

    # ax.legend()

    #scatter plot

    ax.scatter(df['location-lat'], df['location-long'],df['event-id'], c="#1292db", lw=0, alpha=1, zorder=5, s = 3)

    ax.set_xlabel('Latitude')

    ax.set_ylabel('Longitude')

    ax.set_zlabel('Time of Year')

    ax.set_title("Migration of Lesser Black-backed Gulls in 2001",fontsize=20);

    plt.show()

# d = pd.read_csv('migration_original.csv')

#91732 is the id of the bird whose trajectory we are interested in.

#feel free to change the id of the bird

k = df_spacetime.loc[df_spacetime['tag-local-identifier'] == 20527]

​

#normalise the time axis

#use the event-id to determine the time stamp

#Note: event-id is in increasing order between a range of numbers,

#so we can easily normalise for a given bird

#In the next post, I will plot using time-stamps

#and that will make things more clear

​

k['event-id'] = (k['event-id'] - k['event-id'].mean())/k['event-id'].std(ddof=0)